### This script calculates standard errors using the results of the tempdisagg package. 
### This script must be run within the tempdisagg package (where it is sourced). 


n <- length(interpqtrs)

se <- list() 

for (gr in names(scf_data_disagg)) {
  for (l in names(scf_data_disagg[[1]])) {

    u <- scf_data_disagg[[gr]][[l]][["residuals"]]
    
    sigma_u <- sd(u)
    
    # autoregressive parameter (for fernandez, its always 0) 
    a <- scf_data_disagg[[gr]][[l]][["rho"]]
    
    A = matrix(0, nrow=n,ncol=n)
    
    for (i in c(1:(n))) {
      for (j in c(0:(n-1))) {
        if (i+j<=n) {
          A[[i,i+j]] = a^j
          A[[i+j,i]] = a^j
        }
      }
    }
    
    V <- A * (sigma_u^2 / (1-a^2))
    
    scf_expected <- seq(1,n,12)
    
    scf_present_ind <- 1:length(scf_expected)
    
    
    C = matrix(0, nrow = length(scf_expected), ncol = n)
    
    for (row in scf_present_ind) {
      C[row, scf_expected[row]] <- 1
    }
    
    B <- V %*% t(C) %*% solve(C %*% V %*% t(C), tol = 1e-40)
    
    X <- scf_data_disagg[[gr]][[l]][["model"]][1:n,]
    
    tL <- solve(t(X) %*% t(C) %*% solve(C %*% V %*% t(C)) %*% C %*% X, tol = 1e-40) %*% t(X) %*% t(C) %*% solve(C %*% V %*% t(C))
    
    I <- diag(nrow(X))
    
    var <- (I - B %*% C) %*% V + (I - B %*% C) %*% X %*% tL %*% (C %*% V %*% t(C)) %*% t(tL) %*% t(X) %*% t(I - B %*% C)
    
    se[[gr]][[l]] <- sqrt(diag(var)) / sqrt(n)
    

  }
}

for (gr in names(se)) {
  se[[gr]] <- bind_cols(data.frame(date = interpqtrs), se[[gr]])
}

se <- bind_rows(se, .id = "cat") %>% 
  select(date, cat, everything()) %>% 
  mutate(cat = factor(cat, levels = dfa_splits(agg_level = "bus_cycle2")$groupnames)) %>%
  arrange(date, cat) 

# save standard errors of levels 
write_csv(se, file.path(output_directory, "se", dfa_file_names(item = "se", group = group, agg_level = agg_level, method = tempdisagg_method, filetag = filetag, scf_ignore = scf_ignore)))

